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Abstract 

For special model problems, Fourier analysis gives exact convergence rates for the 
two-grid multigrid cycle and, for more general problems, provides estimates of the two- 
grid convergence rates via local mode analysis. A method is presented for obtaining 
multigrid convergence rate estimates for cycles involving more than two grids - using 
essentially the same analysis as for the two-grid cycle. 

For the simple case of the V-cycle used as a fast Laplace solver on the unit square, the 
A>grid convergence rate bounds obtained by this method are sharper than the bounds 
predicted by the variational theory. Both theoretical justification and experimental 
evidence are presented. 
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1. Introduction 


Although the heuristic arguments which indicate that multigrid iterative methods lead to fast ellip- 
tic solvers are easy to understand, they cannot predict convergence rates or guarantee convergence. 
There exist methods which provide a means to calculate theoretical multigrid convergence rate 
bounds for a large class of multigrid solvers for positive definite elliptic problems. Unfortunately, 
the motivation and proofs of these bounds are relatively difficult. For special model problems, a 
simpler calculation, based on Fourier analysis, provides exact convergence rates for the two-grid 
algorithm, but the exact analysis becomes too complicated to carry out for more than two or three 
grids. 

In this paper we present a technique for obtaining convergence estimates for any number of grid 
levels, thereby providing an easy method to verify that multigrid is a fast solver (i.e., only a fixed 
number of iterations will solve the problem). Although the method is strictly valid for only a few 
model problems and for only symmetric U-cycles, it provides a general insight into the details of 
the basic interaction between the coarse grid correction and the relaxation. These new estimates 
of the convergence rates for the A;— grid problem are no more difficult to calculate than the exact 
convergence rates for the two-grid problem. 

In Section 2 we define a simple multigrid U -cycle and derive a formula for the iteration matrix 
associated with a multigrid cycle. A brief introduction to the grid independent convergence rate 
bounds given by the variational theory is given in Section 3. In Section 4 it is shown that bounds 
on the k grid convergence rate can be found in terms of the A: 1 grid convergence rate, using the 

same type of calculation as for obtaining the exact two grid convergence rates. Moreover, since the 
two grid convergence rate can be found exactly, we obtain convergence rate bounds recursively for 
any number of grid levels. It is shown that these bounds are sharper than those obtained by the 
variational theory. In Section 5 we show the details of how the results of Section 4 can be used to 
obtain convergence rates for a model problem. Finally, in Section 6, we compare the convergence 
rate bounds predicted by our theory and the convergence rate bounds given by the variational 
theory. The exact three grid convergence rates are compared to our three grid convergence rate 
bounds. 
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2. Notation 


2.1 Multigrid Cycle 

A sequence of uniform grids is given with mesh sizes h k ( k = 1,2,3,...), where h k +i = h k f 2. 
Consider the discrete equations on the h k grid of the form 

A k U k - F k (1) 

where A k is symmetric, positive definite. Given uj^, an approximate solution to (1), the multigrid 
cycle MG for producing an improved approximation uj^, 

( 2 ) 

is recursively defined as follows: 

^ ^ 1 solve (1) by any direct or iterative method, yielding the final result ti^. Otherwise do 

(A) through (E). 

(A) Perform ri relaxation sweeps on (1), resulting in a new approximation u*. 

(B) Transfer (“restrict”) residual from grid h k to grid h k -\. 

^*-1 — “ A k u k ). (3) 

(C) Starting with = 0, update coarse grid approximation. 

“l-’l - MG(i - - “*.«.))• (4) 

(D) Calculate u* — Ufc + /*_ . where /*_, is a suitable interpolation (“prolongation”) from 

grid h k _! to grid h k . 

(E) Perform r 2 relaxation sweeps on (1), starting with u k and yielding the final result ul 1 *. 
This cycle is called V— cycle or 

Note: The relaxations done in steps (A) and (D) need not be the same. 

2.2 Error Analysis 

As in the analysis of the classical fixed point linear iterative solvers, we are interested in analyzing 
the “iteration matrix” associated with the multigrid process. If the error after the i-th iteration 
can be written in terms of the error before the i-th iteration as: 

V t - „W = B t {U„ - u<"‘>) ( 5 ) 
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then B k is the iteration matrix. We are interested in the convergence factor of the V- cycle, e*, 
given by 

e k = inf{e : \\U k - MG{k,v k , F*)|| < e\\U k - v*|| W* € H k } (6) 

= »n/{e : < e||tnfc|| Vw k eH k } 

= \m 

where || • || is some norm defined on H k , the space of all vectors defined on the grid h k . 

To find an expression for the iteration matrix corresponding to a V— cycle, the steps (A)-(E) 
are rewritten in terms of the errors (£rue solution — current approximation) . Given an initial 
error, ej^ = U k - , the iteration matrix is defined recursively as follows: 

If k = 1, then Bi = 0 (the error is eliminated completely on the coarsest grid). Otherwise 

•f* = P) 

where ej^, the error of the improved approximation ttj^, is given by steps (A')-(E'): 


(A*) After relaxation: 

e„ = (G„) r '4 0) . 

where G k is the iteration matrix of the (pre-) relaxation. 

(B 1 ) After residual transfer: 

U k -1 = AjijjJ -1 A k e k . 

(C*) After computing coarse grid correction, since 4-i = U k - k : 


4-1 = ft-ieCi- 


(£)') After interpolating correction to fine grid: 



* r k (AO) _ (1) 1 

e k - h-l( e k-l e k-l) 

e k -/t 1 (/-B k - 1 )4-r 


( E *) After relaxation: 


e { k ] = (E k y^ 


(i) 


where E k is the iteration matrix of the (post-) relaxation. 
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Therefore for 1 < k < m, 

B„ = (£*)'■(/ - /»_,(/ - ( 8 ) 

and 

Si =0. 

3. Variational Theory 
3.1 Notation 

Consider the Euclidean spaces H k = lR nk , k = 1 , 2, •••, m and two types of inner products on 
these spaces, one type is the usual discrete L 2 inner product in d-dimensions defined as 

(“*> v k)k = ^ X^(ujfc),(v fc ),- 

t 

and the other is the “energy” or “operator” inner product (associated with an operator L k ) defined 
as 

(«*,*>*)* = {u k ,L k v k ) k . (9) 

The norms induced by these inner products are denoted by: 

ll«*IU = 

and 

\\ u k\lk = 

Adjoints relative to the {•, •) inner products are denoted by (■) T and adjoints relative to the (•, •) 
inner products are denoted by (■)*. 

3.2 Bounds Independent of the Number of Levels 

The best available multi-level convergence rate bounds come from the analysis in the “variational” 
setting, i.e., where it is assumed that: 

a. The intergrid transfer operators are related by: 

*t l = til i) T 

and 
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b. the “Galerkin” choice of the coarse grid operators is assumed: 

M-i = 

For the complete theory for general domains and W- as well as K-cycles, see [3]. For positive 
definite symmetric elliptic p.d.e.’s in the unit square, typical V -cycle convergence rate bounds 
which are independent of the number of coarse grids are summarized below. For a more complete 
discussion of the types of convergence proofs, see [6]. 

As in the standard variational multigrid analysis, see e.g., [4], [5] and [1], define two operators, 


S k and T k : 

s k = 

(10) 

and 

& 

H 

i 

(11) 


Then S k and T k are orthogonal projections in the (•,•) inner product. Specifically: 

1 . S k S k = S k , T k T k = T k and T k S k = S k T k - 0. 

2. For any u k , 

u k = T k u k + S k u k 

and since S* k = S k and T k = T k (but Sj # S k and T k ^ T k , in general), we also have 

Mil = p fc u fc ||| + 17^12. (12) 

Note: If u* = S k u k , then u k G ranye(/j_ 1 ) ) i.e. it is representable on the i grid. If 
v k - T k v k , then I k ~ 1 A k v k = 0, i.e. v k G nullapace{I k ~ x A k ) and if v k corresponds to an error, its 
residual cannot be distinguished from the null vector on the grid. 

Let e k be the operator norm convergence factor of the k grid V-cycle, see equation (6). 

Theorem 1 If there exist /3i,/3 2 > 0, independent of k, such that for every u k in H k and for all 
k — 2,3, • • ■, 

ll(G*) r 'u t ||£ + Pi\\T k (G k ) r 'u k \\l < IIKil; (13) 

and 

||( J E7jfc) ra «*||| + p 2 \\T k {E k y*u k \\l < lllu.ll'i (14) 

then for every k = 2, 3, • • •, 

e*<(l + A)" l/, (l + W‘ /! - < 15 ) 
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Proof. See [3]. 

In [3], Theorem 1 is proved using recursive formulas of the form: 


sup ( 
0<<<1 


t + Cy-x(l -Q 

1 +Pt ' 


( 16 ) 


where j > 2. 

Although these formulas appear to give grid dependent bounds, notice that cy = 1/(1 + /9), 
independent of j. Therefore these formulas give the same convergence rate bound for the two grid 
cycle as for an arbitrary number of grids. They cannot give sharper intermediate estimates for a 
specific number of grids. In the next section we show that, under special assumptions, one can do 
better than this. 


4. Multi-Level Analysis 

From equations (8) and (11), the two grid iteration matrix, if 2 , is: 

* 2 = (**)* T 2 (G 2 ) r ‘. (17) 

Under certain assumptions about E 2i /£, if and G 2 , it is feasible to explicitly calculate the 
spectral radius of B 2 and hence find exact asymptotic convergence rates for the two grid algorithm. 
Section 5 contains the details of the calculation of the spectral radius of operators of the slightly 
more general form ( E k ) r *(T k - aS k )(G k ) ri where a is a scalar. Notice that the two grid operator 
is of this form with a = 0. In general, the A;-grid V— cycle iteration matrix is given by 

B k = (E k y*(T k + /*_ x B*_ i A*2 , 1 A*) (G* ) ri (18) 

and is not of this form unless, for instance, Bt-i is a constant multiple of the identity (i.e., all error 
components on the k - 1st grid are uniformly reduced). It is possible, however, to obtain estimates 
of the k grid convergence rates by replacing B fc -i by the constant ||J9t_ 1 ||t_i times the identity. 
For convenience, fix the number of relaxations, r,, and define 

Gj = (GjY 1 - 

We assume the Ay’s are symmetric and that Ay, Gy and Ej commute. 

Theorem 2 If r x = r 2 and for every j — 1,2, •• • ,k, Ej = G*j, then Bj is positive semi-definite 
and 

< 15,11, < isjpv + ib, ) e,i,. ( 19 ) 
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Proof. 

Since = (/y„ 1 ) T , Ty = (Ty)* and Sy = ( Sy)*, each By is symmetric with respect to the 

energy inner product (Bj — (By)*). Each Bj is also positive semi-definite. This is most easily 
seen by induction. B 2 is positive semi-definite since (B 2 W 2 , ^ 2)2 — for every W 2 6 B 2 . 

Assume Bj - 1 is positive semi-definite for some j > 2. From formula (18), 

(BjWj y Wj)j = (GjTjGjWj, wj)j + (Gjlj-iBj-iAj^Ij AjGjWj , u/y)y. 

Since ( Tj ) 2 = Ty for arbitrary Wj G By, we can write 


(Gy TjGjWj , toy)y — (TjGjWj , TjG jWj)j 
and, similarly, since 7y _1 = (/y_i) T , 

The right hand side of the previous equation is non-negative by the inductive assumption and thus 
(Bj w j, Wj)j = (TjGjWj, TjGjWj), + (Bj- l (A]ljfUjGjWj),(Aj}j’- 1 AjGjWj))j- l 


> \mcjWjfj > 0 . 


( 20 ) 


Therefore By is positive semi-definite. 

We can also use equation (20) to obtain the lower bound for the || By ||y ’s in equation (19). 

To prove the upper bound in equation (19), let «?y be an arbitrary vector in By. By the definition 
of «B, and Sj, 


“.•“’.If = 


< 


(B,. J Aj‘ l Ij-%G i w J ,AT} 1 lj- l A,S i a, i ) J . 1 


— III Bj-i ||y-i (GjSjGjWj , Wj)j. 


Therefore, 

(BjWj,Wj)j < ((GjTjGj + \\Bj. l \\j_ l G)SjGj)wj,Wj)j. 

This proves the upper bound in equation (19) since, for any symmetric matrix, Dj, 


( 21 ) 


IP, lb- 


sup 

Wj G Hj 

IKII ± o 


(DjWj, Wj)j 
(Wj,Wj)j 
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Let 71 = 0 and define 7* recursively as 


7 fc = |||G;(T fc + 7t-i5 t )G t || t . (22) 

The next theorem asserts that each 7 * is a bound on the fc-grid convergence rate. 

Theorem 3 For every k>2, 

\\B k \h < Ik- 


Proof. 

From the definition of the 7’s, we see that 

III #2 1 2 = 72 

if Py— ill,— 1 < 7y_ 1 , then for each tiy G fly, 

(BjUj , Uj )j < [G j ( Tj + | Bj-\\lj-\Sj)GjUj , uy) y 

< |||T i G J .u i ||^ + 7y-i||5A«ylli 

= (G*(T i + lj - l S j )G j u j ,u j ) j 

and therefore |||fly||y <7y- 

In order to compare these bounds with the theoretical bounds of Theorem 3.1, we define a 
constant, /9, independent of the number of grid levels, such that 


/? = inf 
*>i 


/ 

inf 
y T k & k u k * o 


( 


IKlllMflGmtH 


l ) 


From this definition we see that, for every j > 1, 

l|Gy«y||> + ^lll^yGyUjIly < ||uy||y 

or equivalently, 

I^Gyyyg + 


(23) 


(24) 
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Note: The constant is given for some common relaxation methods and for piecewise (bi)linear 
finite elements (i.e., (bi)linear interpolation for finite differences) in [3]. Equation (23) holds for 
j3 = 0 whenever the iteration matrix of the relaxation, Gy, is a contraction. The convergence rate 
of the V— cycle can be shown to be bounded away from one only when equation (23) holds for a 
positive /3. 

Theorem 4 For all k > 1, 


ik < 1/(1 + p). ( 25 ) 

Proof. 

Equation (25) holds for 71 since 71 = 0 and /? is non-negative. 

Assume 7^-1 < 1/(1 4" P) for 3 2. For any Wj E Hj, equation (24) implies that 

(G* {Tj + , wj), = {G'TjGjWj, wj ) y + 7,-1 (GySyGytuy, ti>y)y 


Since G* (Tj + 7 y_ 1 S J )G J - is symmetric, then 7, < 1/(1 + P) by definition (22). 

Finally, we note that this new technique has another useful feature. If the fc-grid convergence 
rate is known, then the k + 1-grid rate can be estimated. Thus, if the convergence rate of a V-cycle 
for a certain number of grids is known, it is possible to predict the effect of adding one additional 
grid. 

In the next section we show how to calculate the convergence rate bounds, 7*. 

5. Details of the Fourier Analysis 

For certain model problems, it is possible to calculate exact convergence rates of the two grid 
V-cycle. The same techniques can be used to compute the fc-grid convergence rate bounds of 
Theorem 4 for the symmetric V-cycle. 

The Fourier mode ezp(t£ • i/h*) on level k appears as the mode exp(i’2£ • x//ijt_i) on level k— 1 
since h k = ,5h fc _ 1 . Therefore, on level fc - 1 it coincides with every Fourier mode exp(i9 [ • xjh k ) for 
which 0 ! = 6{modn). Thus restriction operators introduce coupling between each lower mode 0 and 
its (2 d - 1 ) high-frequency harmonics {£' : jt/2 < |£| < *,£! = 9{modv)}. Interpolation introduces 
coupling among the same modes. This is summarized in the following two formulas 


I k -iexp(i9_-xjh) = Y, ^-i(£') ex P(*!' 2/*) ( 26 ) 

0^=0(m.odir) 
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and 


l exp{%0^ • xjh) = 1 (0 , )exp(i0-z/h). 0!_=0(modi r). (27) 

/jt _1 an< ^ ^fc-i(£) are ca H e d the symbols of Z* -1 and i^” 1 respectively. 

Since the fine and coarse grid operators, as well as the relaxation, do not introduce coupling of 
more Fourier modes we conclude that the two level cycle has the 2 d harmonic modes as its invariant 
subspace. Hence the convergence properties of the multigrid cycle can be studied by looking only 
at the interactions of each set of 2 d Fourier components. 

Interpolation and restriction can be represented in matrix form as 

' iLid 1 ) ' 

(£)= : (28) 

. 1(£ 2 ) . 2 J xl 

and 

'f“‘®')] lxJJ - (29) 

The fine grid operator is represented as 

Mi) = 

Mi 2d ) 

the coarse grid operator as 

A k -i(2£)=[A t _ 1 (2£ 1 )] ixi , (31) 

and the relaxation matrix as 

‘ 9(l l ) 

G*(fi)= ••• ^ - (32) 

^(— ) - 2 d x2 d 

In terms of these matrices the two level cycle is represented by the following 2 d x 2 d matrix 

BiW = [O*®]' 2 [/ - ?»-i(fl(Ai-,(2£|)- 1 /*->(«)A 1 (0] [G*®]' 1 • (33) 

Take the simple case of the V'-cycle used as a fast Laplace solver on the unit square. A few 
simple relaxation methods (such as Richardson or damped Jacobi) can be easily analyzed for this 
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particular problem simply because the Fourier modes are eigenvectors of the relaxation. For this 
2 dimensional problem, the two grid convergence rates are obtained by computing the largest (in 
absolute value) eigenvalue over all the 4 X 4 matrices given by formula (33). The bounds for the 
convergence rates for more than two grids are computed in exactly the same manner, inserting a 
constant in formula (33), preceding the interpolation symbol, in order to obtain the 7* of Section 4. 

For completeness we include some examples of some operators and the elements in their symbols. 
If — ^ corresponds to the five point Laplacian, given by the stencil 


A*) - 1 

* 


-1 

-1 4 -1 

-1 




then 


X^(9) _ 4(sin 2 (g 1 /2) + sin 2 (fl 2 /2)) 

k £5 1 


(34) 


(35) 


where £= (0i,0 2 ). 

For bi-linear interpolation and the residual restriction operator which is its adjoint (usually 
called full weighting), the elements of the symbol matrices are given by 




and 


n-'w = (i±^. )( i±^). 


(36) 

2 M 2 ( 37 ) 

A simple relaxation, damped Richardson (which is the same as damped Jacobi in simple cases 
we consider), the elements of G*(£) are given by 

0(£) = 1 - 2 wcA*(£) (38) 

where c _1 = sup*A t (0) and 0 < w < 1. 

When the coarse grid operators are determined “variationally” , i.e., A k -i = l£~ 1 A t l£_ 1 , then 

4 

k * L . * • At „ « 

(39) 


Vi(2£) = Eti(r')A l (tf i )/*- 1 (£ i ). 
1=1 


h 4°' is a nine point discrete Laplacian, given by the stencil 


A m = _L 

* 3 H 


-1 

-1 

-1 

-1 

00 

-1 

-1 

-1 

-1 
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(40) 


hk 


then 


(41) 


« (9 ) 8 - 2 008(0!) - 2cos(0 2 ) - 4cos(0i) cos(0 2 ) 

A k ( 2 )- 3^2 


If is given by the corresponding stencil on the ht-i grid: 


A W - 1 


" 8 * 2-1 


-1 -1 -1 

-1 8 -1 

-1 -1 -1 




(42) 


then 

. f9 ) 8-2 cos(20i) - 2 cos(2fl 2 ) - 4003(20!) cos(2g 2 ) 

This coarse grid operator, though corresponding to the same stencil as the fine grid operator, is in 
fact also the “variational” coarse grid operator when using bi- linear interpolation and full weighting. 
This is convenient when trying to use the variational theory but trying to keep the symbols of the 
discrete operators simple. 

6. Results and Comparisons 

In [2], it is observed that the variationally derived bounds are too pessimistic, at least for typical 
parameters in a multigrid V-cycle when using a modest number of grids. Although in Section 4 
we prove that the fc-grid bounds are smaller than the grid independent bounds of the variational 
theory, it is not known, in general, how much better the new bounds are. The limit of the bounds, 
1 k , ask goes to infinity, does not always degrade to the theoretical rate of 1/(1 + P). For a modest 
number of grids, the new estimates are much closer to the actual convergence rates. 

Moreover, we note that this new technique has another useful feature. If the fc-grid convergence 
rate is known, then the k + 1-grid rate can be estimated. Thus it is easy to predict the effect of 
adding an extra grid to a given multigrid cycle. 

Consider the simplest case of a one dimensional problem given by 

= / on (0, 1) (44) 

u(l) = 6 


_ fu 

dx 2 

u( 0) = a ; 
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Figure 1: The convergence rate bounds for the one dimensional model problem 


Define a V^-cycle using the standard three point discretization of the second derivative, linear 
interpolation and damped Jacobi (with uj = 2/3) for the relaxation. The convergence rate bounds 
predicted in Section 4 for a V(l , 1) cycle (one pre- and one post-relaxation) are given in Figure 1. 
The 7 * are shown as a function of the number of grid levels, k. In contrast, the variational theory 
gives a convergence rate bound equal to 0.408, see [2]. 
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Tables 1 and 2 compare the exact three grid convergence rates (computed with a three grid 
Fourier analysis) to the bounds predicted by our theory. For the optimal w, u = 2/3, the bounds 
stay close for r\ = = 1, • • • ,4. 


ri 

r 2 

u> = 1/2 

u = 2/3 

u = 3/4 

i 

i 

.2665 

.1655 

.2499 

2 

2 

.1086 

.0826 

.0903 

3 

3 

.0743 

.0562 

.0539 

4 

4 

.0570 

.0430 

.0392 


Table 1: One-dimensional exact three grid convergence rates of V(ri,r 2 ) using damped Jacobi with 
damping parameter u 


r\ 

Tl 

uj — 1/2 

u = 2/3 

u = 3/4 

i 

l 

.2667 

.1667 

.2499 

2 

2 

.1161 

.0887 

.0984 

3 

3 

.0804 

.0610 

.0582 

4 

4 

.0617 

.0466 

.0425 


Table 2: One-dimensional estimated three grid convergence rates of using damped Jacobi 

with damping parameter u> 
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Finally, for Poisson’s equation with Dirichlet boundary conditions on the unit square, we com- 
pare our bounds to the asymptotic convergence rates seen experimentally. Using the grid sizes 
indicated in the first column of table 3, we ran experiments using a damped Jacobi relaxation, the 
nine-point discretization of the Laplacian, as given in Section 5 and bilinear interpolation and 
full weighting. Starting with a random initial error, the V^l, 1) multigrid cycle was used, rescaling 
the error after every cycle in order to see the asymptotic convergence rate. The center column 
contains the bounds given by our method. The grid independent bound given by the variational 
theory is 0.40, see [3], 


grid sizes 

our 

bounds 

worst case- 
experimental 

h =(1/4, 1/2) 

.110 

.110 

h = (1/8, 1/4, 1/2) 

.217 

.211 

h = (1/16,1/8,1/4,1/2) 

.258 

.241 

* = (1/32, 1/16, 1/8, 1/ 4, 1/2) 

.286 

.246 

h = (1/8, 1/4) 

.206 

.206 

*=(1/16, 1/8, 1/4) 

.254 

.239 

* = (1/32,1/16,1/8, 1/4) 

.284 

.245 

h = (1/16, 1/8) 

.238 

.238 

h= (1/32,1/16,1/8) 

.275 

.244 


Table 3: Comparison of bounds with actual rates for two-dimensional Laplacian (9 pt. stencil) 
using damped Jacobi with damping parameter .75 and r = 1. The grid independent bound is 0.40. 
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